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A Lattice Boltzmann Relaxation Scheme for Inviscid 

Compressible Flows 

S.V. Raghurama Rao* Rohan Deshmukh^ and Sonrabh Kotnala^ 

Department of Aerospace Engineering, Indian Institute of Science, Bangalore-560012, India 

A novel Lattice Boltzmann Method applicable to compressible fluid flows is developed. This 
method is based on replacing the governing equations by a relaxation system and the 
interpretation of the diagonal form of the relaxation system as a discrete velocity Boltzmann 
system. As a result of this interpretation, the local equilibrium distribution functions 
are simple algebraic functions of the conserved variables and the fluxes, without the low 
Mach number expansion present in the equilibrium distribution of the traditional Lattice 
Boltzmann Method (LBM). This new Lattice Boltzmann Relaxation Scheme (LBRS) thus 
overcomes the low Mach number limitation and can successfully simulate compressible 
flows. While doing so, our algorithm retains all the distinctive features of the traditional 
LBM. Numerical simulations carried out for inviscid flows in one and two dimensions show 
that the method can simulate the features of compressible flows like shock waves and 
expansion waves. 


I. Introduction 

The lattice Boltzmann method (LBM) has emerged as an alternative to traditional CFD methods in the 
recent years. Its simplicity, ability to simulate complex fluid flow phenomena and amenability to parallel 
programming are among the reasons for its rising popularity. The lattice Boltzmann method is based on 
microscopic models and mesoscopic kinetic equations.® The kinetic equations are simplified in the sense that 
the velocity space is discrete instead of the continuous velocity space used in kinetic theory. The lattice Boltz¬ 
mann model thus consists of a small number of fictitious particles streaming (along certain fixed directions) 
and colliding on a lattice. From this simplistic model, the macroscopic equations of fluid flows can be recov¬ 
ered by Chapman-Enskog expansion and thus such a simplistic model can be used to simulate fluid flows. 
This connection between the mesoscopic and macroscopic levels allows LBM to simulate continuum fluid 
flows. A key feature in favour of LBM is that the convection term of the lattice Boltzmann equation is linear 
as opposed to the nonlinear convection term in the macroscopic equations. Extensive reviews of the lattice 
Boltzmann method can be found in the papers of Chen and DoolerP as well as that of Benzi et aP and in 
the books of Rothman & Zaleski,®Wolf-Gladrow,® Rivet & Boon,® Succi,® Sukop & ThornPand Mohamad.® 

Historically, the lattice Boltzmann method was developed as an advancement over the lattice gas automata 
methodiHlEniEI] Later, He and LucJ^^ showed that the lattice Boltzmann equation can be derived directly 
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from the Boltzmann equation. Irrespective of its numerous advantages, the lattice Boltzmann method has a 
serious drawback in that it is restricted to simulation of low Mach number (essentially incompressible) flows. 
There have been many attempts to extend this method to simulate compressible fluid flows. Alexander et 
have attempted to simulate Burgers equation with shocks using LBM. They have modified a parameter 
in the equilibrium distribution so that sound speed can be chosen arbitrarily. A major limitation of their 
model is that it can only simulate an isothermal fluid.^ Shouxin et al^ have used a 13-bit model based on 
the FHP lattice with the particles moving along the links being divided into two types, each type having its 
own energy level and the rest particle having a different (third) energy level. The equilibrium distribution 
is required to satisfy flux conditions besides the conservation of mass, momentum and energy. They have 
simulated the shock tube problem. Guangwu et have used a similar approach on a square lattice and 
have tested their model on three compressible flow test cases. The drawback faced by both the models 
is that too many parameters have to be chosen. SurP^ introduced a semi-discrete LB model wherein the 
particle velocities are determined by the mean flow velocity and internal energy. The particle velocities are 
adaptive which permits the mean flow to have a higher Mach number. The drawback of this model, as 
reported by Sun, is that the relaxation parameter has to be set equal to one. Kataoka and TsutahareESEU 
developed a lattice Boltzmann model for the compressible Euler and Navier-Stokes equations wherein the 
specific heat ratio can be chosen freely. Their model can simulate compressible fluid flows well but is unsta¬ 
ble for Mach numbers greater than unity. Tolk^i^ presented a new model based on the lattice Boltzmann 
method for simulation of thermal compressible flows. He used the lattice Boltzmann equation to solve for 
the flow field while the temperature equation is solved by a finite difference scheme. The sound speed varies 
proportionally with temperature to allow large density variations. The drawback of this model is that it is 
restricted to low Mach numbers while the thermal flows with large density variations can be simulated. Yan 
et alW^ proposed a Lagrangian LBM to simulate compressible isentropic flows. They have used displace¬ 
ment distribution functions instead of velocity distribution function. While they have successfully simulated 
compressible inviscid flows, the authors report that issues regarding accuracy, stability and wall boundary 
condition for their method need to be solved. Zhang et and Yan et have proposed a three-energy 
level and three-speed Lattice Boltzmann model using higher moments of the equilibrium distribution to sim¬ 
ulate compressible flows. Other author^MHHBlIl have used finite difference LBM with a modified equilibrium 
distribution function to overcome the low Mach number restriction of the traditional LB method. Qu et 
have developed an alternative method to construct the equilibrium distribution function to overcome 
the low Mach number restriction of the traditional LBM. 

Although the list presented above is not exhaustive, it can be seen that attempts have been made in several 
directions to develop a robust LB method for compressible flows and yet, a widely-accepted route to tackle 
compressible flows using the lattice Boltzmann framework has not yet emerged. In this paper, we have 
attempted to construct a compressible lattice Boltzmann method using a hitherto unexplored area, namely 
the framework of relaxation systems for hyperbolic conservation laws. This Lattice Boltzmann relaxation 
scheme (LBRS) is based on the interpretation of the diagonal form of the relaxation system as a Lattice 
Boltzmann equation. As a result, the local equilibrium distribution in our method is an algebraic function 
of the conserved variables and the fluxes. LBRS simulates compressible flows successfully as the local 
equilibrium distribution function is not based on the low Mach number expansion. The following section 
formulates the lattice Boltzmann relaxation scheme. In section III, compressible flow simulations for one 
and two dimensions are carried out using LBRS. Section IV concludes this paper. 
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II. Lattice Boltzmann Relaxation Scheme (LBRS) 

The lattice Boltzmann relaxation scheme is based on the framework of a relaxation system. In this section, 
we first utilise the relaxation system of Jin & XirP^ for a hyperbolic conservation law in one dimension 
to formulate our method. This relaxation system is then diagonalized and further interpreted as a discrete 
velocity Boltzmann equation, following Natalini.l^ Then, we extend LBRS to two-dimensions. The relaxation 
system of Jin and Xin cannot be diagonalized for the case of two-dimensions. We use a novel diagonal 
relaxation system which is isotropic to extend our scheme to two-dimensions. 


A. LBRS for scalar equations 

The relaxation system was introduced by Jin and XirP^to replace a nonlinear conservation law with a semi- 
linear hyperbolic system. We demonstrate the interpretation of this relaxation system as a lattice Boltzmann 
equation first for a one-dimensional hyperbolic conservation equation like the inviscid Burgers equation. 
Consider the inviscid Burgers equation. 


^ , dgju) 
dt dx 


= 0 , 


{x,t) G X R+, u G R^ 


( 1 ) 


with initial data u{x, 0) = Uq{x) and g{u) = 

The relaxation system introduced by Jin and XiiP^ is given as: 

du dv 


dt dx 


= 0 , 


V G R^ 


dv ,odu 1, , ,, 


( 2 ) 


where e is the relaxation parameter and A is a positive constant to be determined from the sub-characteristic 
condition. Note that as e —> 0, the second equation of the relaxation system gives v = g{u) which, when 
substituted into the first equation of the relaxation system leads to the original conservation equation 0 . 
The relaxation system can be written in vector form as follows: 




( 3 ) 




U 


0 1 


0 

where 

Q = 

V 

; A = 

1 

> 

to 

O 

1_ 

; H = 

-^{v-g{u)) 


Note that the eigenvalues of A are real and distinct (±A) and thus the relaxation system is strictly hyperbolic. 
Therefore, the diagonal form of the relaxation system is obtained by introducing characteristic variable vector 
f = i?“^Q in 0 where R is the right eigenvector of A, given by the following expression. 
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The vector form can be written in terms of the components as 


dt dx 



5(w)] 


df2 , % ^ 

dt dx 



9{u)] 


(5) 


This system can be interpreted as a discrete velocity Boltzmann equation if we introduce the following new 
variables: 


req 

2 2A 


and 


feq^u g{u) 
" 2 + 2A 


( 6 ) 


Equation now transforms to 


dfi_.dh 

dt dx 
9f2 ^ ydf2 




(7) 


with fi and /2 given by equation]^ Thus, —A & A are the discrete velocities and fi & /2 are the corresponding 
components of the distribution function f in the discrete velocity Boltzmann equation with the B-G-K model 
given by 


This interpretation was given by NataliniP^ The initial condition u{x,0) = uo{x) can be rewritten as 
fi = ff‘^{uo{x)). Equations 0 resemble a two-velocity lattice BGK model with the distribution function f 
relaxing to the equilibrium distribution function f®'* within the relaxation time e. This interpretation of the 
relaxation system is the basis of the Lattice Boltzmann Relaxation Scheme (LBRS). In accordance with the 
terminology adopted for lattices in the standard LB literature, this model can be termed as the D1Q2 model 
with two particles traveling in opposite directions with speed A at every lattice node as shown in figure 0 . 
It is to be noted that the equilibrium distribution functions in LBRS defined by equation 0 are simple 
algebraic functions of the conserved variable and the associated flux. These functions are not expressed 
as a polynomial (Gaussian) function of the macroscopic velocity (as is the case in the traditional lattice 
Boltzmann method) and hence our method is able to simulate compressible flows as no low Mach number 
expression is involved. The left hand side of the equation represents convection on the lattice while the right 
hand side represents collision at a lattice node. The LBRS algorithm retains the hopping-collision cycle of the 
traditional method. The variables u, v and g{u) of the relaxation system are recovered as a simple summation 
of the distribution functions as illustrated in equation 0. The original conservation law is then recovered 
from the relaxation system as the relaxation time tends to zero. Thus our method successfully overcomes 
the low Mach number restriction plaguing the traditional LB method while maintaining the structure and 
ease of programming of the LBM. 


u = /i+/2 = /r + /2^ ^^ = -A/i+A/ 2; 5(n) =-A/r« + A/2"« (9) 

Another feature of the traditional LB method that carries over to LBRS is that the number of particles in 
the model can be increased as long as the moment relations are satisfied. Generalizing to n particles, the 
moment relations are 

n n n n 

= = = AJf'^ (10) 

i—1 i—1 2=1 2=1 
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Figure 1. D1Q2 model 


The speed of the lattice particle (A) is determined from the sub-characteristic condition defined in Jin and 
XinP 


> 


dgju) 

du 


or 


^ dg(u) , 
- A < - 4 ^ < A 
ou 


( 11 ) 


The equation Q can be written in discrete from as follows: 

fi{x + Aa;,t + At) = - — [f^{x,t) - f4{x,t)\ ; fori = 1,... 


( 12 ) 


The time interval At is chosen such that convection is exact ie., in the time interval At, the particles 
hop exactly from one lattice node to the adjacent node in the corresponding direction. The discrete form 
(equation fill) can thus be written as 


fi{x + AjAt,t + At) = (1 - uj)f^{x,t) + ujf^‘^{x,t)] ; fori = 1, 


(13) 


where oj = At/e and takes a value between 0 and 2. The limits for the value of uj are decided based on 
a Chapman-Enskog type analysis of the LBRS equations. We present the analysis for the case of one¬ 
dimensional Burgers equation in the following subsection. 

B. Chapman-Enskog analysis of LBRS 

The macroscopic equations are derived from the traditional lattice Boltzmann equations by carrying out a 
multi-scale analysis which is similar to the Chapman-Enskog analysis in deriving Navier-Stokes equations 
from the classical Boltzmann equation. As a result of this Chapman-Enskog type analysis, the diffusive term 
and the corresponding transport coefficient can be recovered. Since the relaxation approximation is known 
to be a vanishing diffusion model to the original governing equation, we have carried out a similar multi¬ 
scale analysis for our scheme, for the case of one-dimensional Burgers equation, to determine the diffusion 
coefficient. The limits for the relaxation parameter are then decided based on the physically apt values of 
the diffusion coefhcient. 


Consider the discrete LBRS equation given by equation (12). 


Aa;,t-b At) = /,(a;,t)- ^[f^{x,t) - f4{x,t)\ ; 


for i = 1,..., n 


(14) 


Substituting Taylor series expansion for the term on the LHS of this equation and imposing the exact hopping 
condition (A^At = Ax) results in 


. 9 f. , 9 f. , ^fAt^ , At2 av. , , ^,2 9y^ , , n 


(15) 


Introducing the time and space scales corresponding to relaxation (tr, x), convection (ti = e Xi = e ^cc) 
and diffusion (t 2 = e“^tr-, X 2 = e~^x) phenomena, the temporal and spatial derivatives are given by 


d 


d 


d 


dt ^ dti ^ dt2 


d d 
dx dxi 


(16) 
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+ 2e^ 


dtidt2 


+ £^ 


dt‘^2 


2 3 2 9^ 

dtdx ^ dtidxi dt2dxi ’ dx^ ^ dx'^i 

Based on the same scalings, we expand fi as 


/. = /r + 


■ e/f ^ + eVi 


( 2 ) 


+ 0{e^) 


(18) 


(19) 


Upon substituting the scaled time and space derivatives as well as the expansion of the distribution function, 
we obtain the following expressions 


0(e): 


5 /r ^ 


dti 


dxi 


il/W=0 


( 20 ) 


O(e^) ■■ 


dfl 


( 1 ) 


+ ^ + A, 


(1) 




At 92/r 




( 21 ) 


9ti ' 9^2 '* 9xi 2 dx'^i 9xi9ti 2 9t23 

We have separated the 0(e) and the 0{£^) terms and hence we obtain two expressions. We now take 


moments of these terms based on the moment relations given by equations (10). Noting that the moments 


of both the distribution function and its equilibrium part must be the same, the zeroth moment of the 0(e) 
terms results in the following expression 



= 0 


( 22 ) 


The zeroth moment of the 0(e2) terms gives 


r2^=e2At| 
dto 




92 

V- V 

9a;2i 

dx'^i 


(23) 


Upon addition of equations (22) and (23) and using the definitions for the scaled time and space derivatives 


we recover the original macroscopic equation with a diffusive correction as shown below. 


du dg(u) . / 1 1 


i=l _ 

9a; 2 


92 

9a;2 


(24) 


Here g{u) = We can conclude from the previous equation that uj must take values in the interval (0, 2) 
to ensure a positive viscosity coefficient. 


C. LBRS for one-dimensional system of equations 


Now that we have discussed our method in detail for the one-dimensional scalar equations, we extend it to 
systems of conservation laws in one and two dimensions. In this subsection, we demonstrate LBRS for 1-D 
Euler equations which are given by 


9U 9G(U) 
dt dx 


( 25 ) 


































where 


U = 


P 

pu 

pE 


is the vector of conserved variables and 


G(U) = 


pu 

p + pu^ 
pu + puE 


P 1 

is the flux vector. E = —-^ H—is the total energy and the equation of state is given by p = pRT. 

^ . . . . . r-i 

The equations for the lattice Boltzmann relaxation scheme for this case are similar to equations Q except 

that the distribution functions are now vectors themselves. This is because each of the mass, momentum 
and energy conservation equations is now represented by a 2-velocity diagonal relaxation system. The LBRS 
equations for the D1Q2 model for the one-dimensional Euler system are 


5fi , 5fi _ 1 

dh .dh _ Ir. .eq. 


(26) 
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where fi = 

/l2 

; ^2 — 
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J 22 
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The moment relations are given by 
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n 

peg 

Jil 
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u = 

pu 

= E 
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= E 

peg 
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; G(U) = 

P + pu^ 

= E^* 
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. fiS . 
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where n is the number of particles. The expressions for the equilibrium distribution functions are 


feg _ 

Jl.k ~ 


Ufc G(U)k 


2 2A ’ 2 ' 2A 


_ 

J2.k ~ 


Ufc , G(U)k 


k = 1,2,3 


(27) 


(28) 


The speed A of the particles is calculated from the sub-characteristic condition specified by Jin and XirPS 
for systems of conservation laws. Out of the two choices proposed by them, we use the following one: 


A = TOoa :[|m — a|, |w|, |m-I- a|] (29) 

We demonstrate in the next section that the D1Q2 (2-velocity) model, though it simulates the features of 
1-D Euler test cases, needs some improvement. Thus we introduce here the D1Q3 model by adding a rest 
particle to the 2-velocity model. We then test this model for the one-dimensional Euler case. The D1Q3 
model, (figure which is the standard model for 1-D flows in the traditional LB method, overcomes the 
flaws exhibited by the D1Q2 model as will be elaborated in the next section. 
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Figure 2. D1Q3 model 


The equations for the D1Q3 model are: 


afi > 5fi _ 1 

9^2 _ 1 rn neqi 

5f3 Sfa _ 1 eq. 

For simplicity, we write the above equations in diagonal form: 


( 30 ) 


(31) 


where f and are 3-dimensional vectors and A is a 3 x 3 diagonal matrix given by 
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1_ 
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1 
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The moment relations for this model remain unchanged from equations (27) except that n (number of 
particles) is now equal to 3. The unchanged moment relations allow us to recover the governing equations 
(in the zero relaxation limit) irrespective of the number of particles used in the model. Addition of an extra 
particle however alters the equilibrium distribution functions. The equilibrium distribution functions for the 
D1Q3 model are evaluated from the following moment relations: 




G(U), = ^A/);^ 


A:= 1,2,3 


(32) 


Assuming the equilibrium distribution functions to be a linear combination of Uj, and G(U)k we obtain 
expressions for fc = 1, 2, 3, subject to conditions (32) 


req _ Ufc _ 

/I.fc - 3 


G(U) 


k . 


req Ufc 

j2,k 3 ’ 


req _ rc 

2A ’ 3 ’ “ T ' 2A ’ 

We now move on to demonstrating the extension of LBRS to two dimensions. 


Ufc , G(U), 


k = 1,2,3 


(33) 


D. LBRS for t-wo-dimensions 


As an improvement over the relaxation system introduced by Jin and XirP^ for two-dimensions, Aregba- 
Driollet & NatalinP^ introduced diagonal discrete velocity systems in 2-D. However, these systems are not 
isotropic and we use the isotropic (diagonal) relaxation system introduced by Raghurama Rao and used by 
Jayaraj^and Arun et This is a D2Q4 model with four particles travelling in the diagonal directions 

with speeds equal to v^A at each node (see figure ^). Let us consider a two-dimensional scalar conservation 
law: 


^ dgiju) dg 2 (u) ^ 

dt dx dx 


(34) 
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The equations for the D2Q4 model which is a relaxation approximation for this scalar conservation law are 
as follows: 


dt dx dy 

dh _ ' ^/2 

dt dx dy 

dh , X ^/s , i dh 
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dt dx dy 
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Equations 1^ can be written in vector form as 
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L /4 J 


The equilibrium distribution functions are determined from the moment conditions: 


(35) 


(36) 


M = Pf = PP‘i; 5i(w) =PAiP'i; g2{u)=PK^i^^ 


(37) 


where P = [1 1 1 1], u is the conserved variable, gi(u) and g 2 iu) are fluxes in the x and y directions 
respectively. 



f 


3 


f 
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Figure 3. D2Q4 (Isotropic relaxation system) model 


Again assuming the equilibrium distribution functions to be a linear combination of the conserved variable 


and the fluxes, we obtain expressions for the equilibrium distribution functions subject to conditions (37). 
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To recover the relaxation system corresponding to the D2Q4 model (the isotropic relaxation system), we 


multiply equations [35] by P, PAi and PA 2 respectively to get 

dv2 


du dvi 
dt dx 


dy 


= 0 


dvi , o du dw 1. , ,, 

dv 2 dw ,odu 1, , ,, 


(39) 


where vi = PAif, V 2 = PA 2 ^ and w = PAiA 2 f = PA 2 Aif = 0. Thus, we recover the relaxation system of 
Jin and XinP^ In the limit e 


0 the original conservation law (equation (34)) is recovered. 


The diagonal relaxation system for a system of conservation laws is similar to equation (35) except that the 


components of the distribution function vector are vectors themselves be., each equation of the system is 
represented by a diagonal relaxation approximation given by equations ( [3^ . For a multi-dimensional system 
of hyperbolic conservation laws, the stability condition proposed by BouchuiP^ must be used. This condition 
states that the eigenspectrum of M (u) must be positive, be., 


cr^M (u)^ C [0,-|-oo[ 


(40) 


If condition (40) is satisfied then there exists a kinetic entropy for the relaxation approximation associated 
with any convex entropy ry of the original conservation law and Lax entropy inequalities are satisfied in the 
limit e —)■ 0. Condition (40) for the D2Q4 model simulating the two-dimensional Euler equations gives 


A = 


max pi 


, -I- u|, |u -I- w -I- ^/l2a)\, \u + v - ^/l2a)\, |u - u|, |u - u - ^/l2a)\, \u - v + A/(2a)|^ 


(41) 


where u, v are the macroscopic velocities in the x, y directions respectively and a is the speed of sound. 


For a typical multi-dimensional discrete velocity Boltzmann equation of the type 


9/. V a/, _ 1 


the LBRS equation is given by 


fi (b) -k AiAt, t + Atj = (1 - w) fi (r), t) + uif^‘ 


where uj = —. The solution, as in the 1-D case, consists of a streaming step and a collision (relaxation) 
step, with streaming in the appropriate directions dictated by the discrete velocity model (in our D2Q4 case, 
the diagonal directions). 


E. Implementation of wall boundary conditions 

For the lattice Boltzmann relaxation scheme presented in the previous section, we need to add a strategy for 
implementation of the wall boundary conditions for our method. We will discuss these boundary conditions 
for the inviscid two-dimensional Euler equations. As the flow is inviscid, the appropriate wall boundary 
condition is the free-slip or flow tangency boundary condition. In the traditional lattice Boltzmann method, 
the principle of specular reflection is used at the wall to enforce the free-slip boundary condition.!^ We 
incorporate this principle in our framework. 
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The free-slip boundary condition implies that the normal component of the velocity at the wall is zero and 
that the flow is purely tangential at the wall. In our framework, the normal component of velocity (actually 
the flux involving the normal velocity component) is expressed as a moment of the distribution functions. 
Let us consider an upper wall (figure Q) where the free-slip condition is to be imposed. 


f. 


f 



Figure 4. D2Q4 model at an upper wall 


For the continuity equation we have, pv = PA 2 P = A(—/f — /2 + /s + fi) = 0 
This condition is satisfied when 

ft = ft 


(42) 


ft = ft (43) 

The superscript a is used to identify the distribution functions for the continuity equation. As a result of 


equations 42 and 43 we have determined the unknown distribution functions (fi-, f 2 ) at the wall while satisfy¬ 


ing the free-slip condition. Similar strategy is followed to determine the unknown distributions corresponding 
to the momentum and energy equations as shown below. 

For the cc-momentum equation, puv = PA 2 f'^ = A(—/f — ft + ft + ft) = 0 
This condition is satisfied when 

fi = ft (44) 


rb _ rb 

J 2 — Js 

For the y-momentum equation, puv = PAiP = A(—/f + ft + ft — ft) 
This condition is satisfied when 

ft = -ft 


= 0 


(45) 


(46) 


/ C PC 

2 ~ ~/3 

For the energy equation, we have Pv + pvE = PA 2 f'^ = ^{—ft — /2 + /a + ft) = 0 
This condition is satisfied when 

ft = ft 


fd _ fd 

J2 — J3 


(47) 

(48) 

(49) 
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III. Results and Discussion 


In this section, results obtained by LBRS for several benchmark test cases for inviscid compressible flow are 
presented. We have included results for both hyperbolic scalar conservation laws (Burgers equation) in 1-D 
and 2-D and hyperbolic vector systems (Euler equations or equations of inviscid compressible flows) in 1-D 
and 2-D. 

A. Sod’s shock tube 

In this test case, gas at two different pressures is separated by a diaphragm. Upon rupturing of the di¬ 
aphragm, an unsteady flow consisting of a moving shock, an evolving simple centered expansion fan and a 
moving contact discontinuity is established. The initial conditions for this test case are 


for a; < 0, 


We have divided the domain [-10, 10] into 50 cells. The results are presented at time t = 0.01 seconds. In 
Figure ([^, the density plot obtained by the D1Q2 model is compared with the analytical solution. The 
result for the D1Q2 model exhibits a staircase-like structure which is similar to the result produced by 
Lax-Friedrichs method in CFD. With the addition of a rest particle, the D1Q3 model however gives smooth 
results (see figure Q). Thus the three-velocity model is an improvement over the two-velocity model, though 
both simulate the compressible flow and resolve all the flow features involved. However, we can conclude 
from the Sod’s Shock tube test case results that both the models involve significant numerical diffusion. This 
is expected because the number of equations is typically more than the number of equations in the original 
conservation laws and numerical diffusion gets augmented for each equation. It may be possible to control 
this numerical diffusion in the LBRS framework but this is beyond the scope of the present work and is a 
possibility for future research. 


PL 


1.0 kg/m^ 


PR 


0.125 kg/mf 

UL 

= 

0.0 m/s 

; forx> 0, 

Ur 

= 

0.0 m/s 

. PL . 


. 100000 iV/m^ _ 


. PR . 


. 10000 iV/m^ _ 


B. Steady Shock 

The next test case presents the LBRS solution for a steady shock in a shock tube for the 1-D Euler equations. 
This test case is taken from Zhang and Shu.^ The computational domain [-1, 1] is divided into 400 grid 
points and the shock is located at a: = 0. The LBRS result is shown in figure Q. The initial Mach number to 
the left of the shock is equal to 2 and the post-shock conditions are determined from the Rankine-Hugoniot 
relations. LBRS (with D2Q3 model) captures the steady shock quite well. 
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Figure 5. Density plot for shock tube test case for the D1Q2 model 



Figure 6. 


Comparison of shock tube test case results for the D1Q3 and D1Q2 models 


C. steady contact 

We tested our scheme on the steady contact test cas^^ for the 1-D Euler equations. The initial conditions 
for this test case are pz, = 1-4, ul = 0.0, pl = 1.0 for x < 0 and pn = 1.0, ur = 0.0, pr = 1.0 for x > 0. The 
domain [0, 1] is divided into 100 grid points and the solution is computed at time t = 2. Figure Q presents 
the result with LBRS for this test case. The contact discontinuity is captured reasonably well, though there 
is significant numerical diffusion present. 
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Figure 7. LBRS result for the steady shock test case 



Figure 8. LBRS result for the steady contact discontinuity test case 


D. Two-dimensional scalar conservation law 


We now present results for a two-dimensional scalar conservation law using the D2Q4 model. We discuss 
two test cases from Spekreijse.l^ The governing equation for these two test cases is the 2-D inviscid Burgers 
equation: 

rl'}i fif'jj^/0') ri^i 

(50) 


du 

Ih 


d{u^/2) du 

——— A - =0 

dx dy 


The problem is solved in the domain [0, 1] x [0, 1] on a 64 x 64 grid. The boundary conditions for the first 
case are: M(0,y) = l; «(!,?/) = —1; 0) = 1 — 2a; 

These boundary conditions result in a normal shock originating at x = 0.5, y = 0.5 and a smooth variation 
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resembling an expansion fan below the shock. Figure Q shows that both these features are captured well 
by LBRS. 



Figure 9. density plot for case 1 of the two-dimensional Burgers equation 


The boundary conditions are altered for the second test case such that an oblique shock occurs in 


right corner a shown in figure (10) instead of a normal shock. The conditions at the boundary for 
are: 


the top- 
this case 


u{0, y) = 1.5 ; u(l, y) = —0.5 ; u{x, 0) = 1.5 — 2x 

Again, both the shock and the smooth variation are captured well by LBRS. 



Figure 10. density plot for case 2 of the two-dimensional Burgers equatuion 

After demonstrating that our method works satisfactorily for a scalar two-dimensional conservation law, we 
now present results for the two-dimensional Euler equations. 

E. Oblique shock reflection test case 

This test case deals with the reflection of an oblique shock wave from a solid wall. The boundary conditions at 
the inflow and the top boundary (post-shock boundary conditions) are set such that an oblique shock enters 
the domain from the top-left boundary. This shock is then reflected from the bottom wall. The values of the 
primitive variables at the inflow boundary are p = 1.0, u = 2.9, v = 0.0, P = 1.0/1.4 and the corresponding 
post-shock values at the top boundary are given hy p = 1.69997, u = 2.61934, v = —0.50633, P = 1.52819, 
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calculated from the oblique shock relations from gas dynamics. The free-slip boundary condition discussed 
in the previous section is applied at the solid wall. The density contours obtained using LBRS are plotted 


in figure (11). It can be seen from the figure that LBRS captures both the incident and the reflected shocks 
satisfactorily. 



Figure 11. Density contours for the oblique shock reflection test case on a 240x80 grid 


F. Explosion problem 

We now present results for the explosion problem which can be considered as a two-dimensional extension 
of the Sod’s shock tube problem.l^ The computational domain is a square of dimensions [-1, 1] x [-1, 1] and 
is divided into 400 x 400 cells. The initial conditions are 


p 


’ 1.0 ’ 


P 


’ 0.125 ’ 

u 

_ 

0.0 

for + y^ < 0.4^ ; 

u 

_ 

0.0 

V 


0.0 


V 


0.0 

p . 


. 1-0 . 


. p . 


0-1 


for +y^ > 0.4^ 


We plot the results for time t = 0.25. At this time, the solution consists of a circular shock wave traveling 
away from the origin, a circular contact surface traveling in the same direction as the shock and a circular 
expansion wave traveling inwards towards the origin. LBRS captures all the three waves accurately (with 


respect to position) as shown in figures (12), (13) and (14). 
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Figure 12. Density contours for the explosion case 



Figure 13. Density along the centreline y = 0 for the explosion case 
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Figure 14. Pressure along the centreline y = 0 for the explosion case 


for X >0, y > 0 ; 


for cc < 0, 2 / > 0.0 ; 


G. 2-D Riemann problem 

The two-dimensional Riemann problem is solved on the square [0, 1] x [0, 1]. The domain is divided into 
four quadrants, each of them having different initial conditions. The initial conditions are: 

1.1 
0.0 
0.0 

1.1 

1.1 

0.8939 
0.8939 

1.1 


p 


u 


V 


_ p . 

- 

p 


u 


V 


_ p . 



p 


0.5065 

u 


0.8939 

V 


0.0 

_ p . 


0.35 


for X < 0, y < 0 ; 


P 


0.5065 

u 


0.0 

V 


0.8939 

_ p . 


0.35 


for X > 0, y < 0.0 


The test case presented here corresponds to case 4 in Liska and Wendroff.^^ At each interface, we have 
a one-dimensional Riemann problem. The initial conditions are selected in such a way that a single wave 
(either a shock, contact-slip or a rarefaction wave) is produced for each 1-D Riemann problem. For the given 
initial conditions, the solution consists of four shock waves evolving in the domain, two of which are straight 
shocks while the other two are curved shocks which border the lens-shaped region (see figure ([15])). The 
LBRS result captures these flow features well as can be seen from the figure (15). Liska and Wendroff report 
that the solution must be symmetric about the axis of the lens and our scheme maintains this symmetry. 
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Figure 15. Density contours for the 2-D Riemann problem on a 400 x 400 grid 


H. Supersonic flow over a forward-facing step 

This test case simulates supersonic flow over a forward facing step in a channel. The channel dimensions 
are [0, 3] x [0, 1] and the inflow Mach number is 3. The initial conditions are set to the inflow boundary 
values and the free-slip boundary condition is enforced at the walls. The density contours (see figure (HI)) 
are plotted at time t = 4. Though the incident and reflected shocks, the expansion fan originating the tip of 
the step and the shock-expansion interaction are captured well, the contours for this test case show a lot of 
non-smoothness. We noticed this non-smoothness in the case of vector conservation laws and only for some 
test cases. An improvement to overcome this non-smoothness is suggested in the conclusion section. 



Figure 16. Density contours for Mach 3 flow over a forward-facing step on a 240x80 grid 


I. Shock diffraction 

In this test case, a moving normal shock (Mach number 5.09) encounters a backward-facing step. Initially, 
the normal shock is placed along the tip of the step. The shock wave diffracts as it passes over the step. The 
computational domain for this case is a unit square [0, I] x [0, 1]. The LBRS simulation is shown in figure 
The density contours are plotted at time t = 0.1561. Here too, the contours are non-smooth but the 
flow features are all resolved well. 
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J. Supersonic flow over a compression ramp 


Result for this test case involving supersonic flow (M=2) over a compression ramp is shown in the (figure [T^ . 
Here too, the flow features involving initial and reflected oblique shocks, expansion wave and shock-expansion 
interaction are captured well but the non-smoothness in the contours appear again. 



Figure 17. Density contours for the shock diffraction test case on a 400x400 grid 



Figure 18. Density contours for supersonic flow over a ramp in a channel on a 240x80 grid 


We have presented results for a variety of test cases simulating inviscid compressible flows to showcase the 
capability of LBRS in resolving flow features. Our scheme successfully captures the flow features involved in 
all the test cases. However, the results for some of the cases for two-dimensional flows, for vector conservation 
laws, show non-smoothness in the contours. We guess that as in 1-D where addition of more discrete velocities 
improved the results, increasing discrete velocities may solve the problem of non-smoothness. Numerical 
experiments with a D2Q9 model will be presented in a forthcoming paper,together with several additional 
simulations for compressible flows with curved bodies. 
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IV. Conclusion 


A simple method based on the interpretation of the diagonal form of the relaxation system as a discrete 
velocity Boltzmann equation and further utilization of this relaxation system as a lattice Boltzmann equation 
has been presented in this work. This novel lattice Boltzmann relaxation scheme (LBRS) does not involve 
any low Mach number expansion of the equilibrium distribution functions as in the traditional LBM. The 
equlibrium distribution functions in this method are simple algebraic functions of the conserved variable and 
the fluxes. Thus, the LBRS can successfully simulate compressible fluid flows. We have tested our scheme on 
a variety of standard test cases both for one and two dimensions and for both scalar and vector conservation 
laws, namely for the inviscid Burgers equation and the Euler equations of compressible fluid flows in both 
1-D and 2-D. LBRS successfully captures all the flow features like shock waves, contact discontinuities and 
expansion waves. In the 2-D cases for some of the test cases for vector conservation laws, non-smoothness 
in the contours is observed and an improved model to over come this drawback will be presented in a future 
worlP^ together with more simulations for compressible flows with curved boundaries. 
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